library(Seurat)
library(Matrix)
library(future)
library(future.apply)
library(ggplot2)
library(dplyr)
library(cowplot)
library(grid)

getGenePlots = function(gene, imgalpha = 0, qmax = NULL, topcolor = "#FF0000", lowcolor = "#FFFFFF", pt.alpha = 1){
  g = SpatialFeaturePlot(DATA, gene, stroke = NA, combine = FALSE, image.alpha = imgalpha, alpha = pt.alpha)
  
  if(!is.null(qmax)){
    qval = quantile(do.call("rbind",lapply(g,function(x){x$data}))[,3],qmax)
    g = SpatialFeaturePlot(DATA, gene, stroke = NA, combine = FALSE, image.alpha = imgalpha, max.cutoff = qval, alpha = pt.alpha)
  }
  
  if(imgalpha == 0){
    g = lapply(g, function(x){x + theme(panel.background = element_rect(fill = "#CCCCCC"), panel.grid = element_blank())})
  }
  genemax = max(do.call("rbind",lapply(g,function(x){x$data}))[,3])
  sfill = scale_fill_gradient(limits = c(0,genemax), low = lowcolor, high = topcolor)
  
  return(lapply(g,function(x){x + sfill}))
}

plotGene = function(gene, imgalpha = 0, extratext = NULL,qmax = NULL, topcolor = "#FF0000", lowcolor = "#FFFFFF", pt.alpha = 1){
  g = getGenePlots(gene = gene, imgalpha = imgalpha, qmax = qmax, topcolor = topcolor, lowcolor = lowcolor, pt.alpha = pt.alpha)
  
  return(plot_grid(get_legend(g[[1]]), plot_grid(plotlist = c(lapply(g, function(x){x + NoLegend()})), ncol = 4), textGrob(extratext), ncol = 1, rel_heights = c(0.2,0.9,0.3)))
  
}

plotGeneColoc = function(gene1 , gene2, extratext = NULL,qmax = NULL){
  g1 = SpatialFeaturePlot(DATA, gene1, stroke = NA, combine = FALSE, image.alpha = 0)
  g1 = lapply(g1, function(x){colnames(x$data)[3] = gene1; return(x)})
  
  if(!is.null(qmax)){
    qval = quantile(do.call("rbind",lapply(g1,function(x){x$data}))[,gene1],qmax)
    g1 = SpatialFeaturePlot(DATA, gene1, stroke = NA, combine = FALSE, image.alpha = 0, max.cutoff = qval)
    g1 = lapply(g1, function(x){colnames(x$data)[3] = gene1; return(x)})
  }
  
  g2 = SpatialFeaturePlot(DATA, gene2, stroke = NA, combine = FALSE, image.alpha = 0)
  g2 = lapply(g2, function(x){colnames(x$data)[3] = gene2; return(x)})
  
  if(!is.null(qmax)){
    qval = quantile(do.call("rbind",lapply(g2,function(x){x$data}))[,gene2],qmax)
    g2 = SpatialFeaturePlot(DATA, gene2, stroke = NA, combine = FALSE, image.alpha = 0, max.cutoff = qval)
    g2 = lapply(g2, function(x){colnames(x$data)[3] = gene2; return(x)})
  }
  
  colors = lapply(1:length(g1),function(x){rgb(g1[[x]]$data[,gene1]/max(g1[[x]]$data[,gene1]),g2[[x]]$data[,gene2]/max(g2[[x]]$data[,gene2]),0)})
  
  g = lapply(1:length(g1), function(x){
    #g1[[x]]$data[,3] = as.factor(rownames(g1[[x]]$data))
    g1[[x]]$data$cell = as.factor(rownames(g1[[x]]$data))
    names(colors[[x]]) = rownames(g1[[x]]$data)
    g1[[x]]$mapping$fill = g1[[x]]$data$cell
    g1[[x]] + scale_fill_manual(values = colors[[x]]) + NoLegend()
  })
  
  
  l1 = get_legend(ggplot(do.call("rbind",lapply(g1,function(x){x$data})),aes_string(fill = gene1)) +
                    geom_point(x = 1, y = 2) + 
               scale_fill_gradient(low = "black",high="red") + theme(legend.position = "bottom"))
  l2 = get_legend(ggplot(do.call("rbind",lapply(g2,function(x){x$data})),aes_string(fill = gene2)) + 
                    geom_point(x = 1, y = 2) + 
               scale_fill_gradient(low = "black",high="green") + theme(legend.position = "bottom"))
  
  return(plot_grid(plot_grid(l1,l2), 
                   plot_grid(plotlist = c(lapply(g, function(x){x + NoLegend()})), ncol = 4),
                   ncol = 1, rel_heights = c(0.2,0.9,0.3)))
  
}

1 Load data

DATA = readRDS("Visium_dataset_irradiation.rds")

DATA$orig.ident = factor(DATA$orig.ident, 
                         levels = c("STD_irrad0", "GW_irrad0",
                                    "STD_irrad3", "GW_irrad3"))

2 H&E images

To reduce file size, the H&E images have been removed from most plots. The images are shown here and can easily be combined with the spatial plots in e.g. Affinity.

g = SpatialFeaturePlot(DATA, "Abca1", stroke = NA, combine = FALSE, alpha = 0)

print(plot_grid(NULL,NULL, plot_grid(plotlist = c(lapply(g, function(x){x + NoLegend()})), ncol = 4), ggplot(x = 1, y = 1) , NULL, NULL, nrow = 3, ncol = 2, rel_widths = c(3,1), rel_heights = c(0.2,0.9,0.3)))

3 cNMF 3 factors

cNMF was performed using the Python package https://github.com/dylkot/cNMF. The results can be found at in the cNMF directory in the same directory as this report.

usages = read.table("cNMF/all_cNMF.usages.k_3.dt_0_5.consensus.txt", header = TRUE)

spectra = read.table("cNMF/all_cNMF.gene_spectra_score.k_3.dt_0_5.txt", header = TRUE)



usages_1 = apply(usages,1, function(x){x/sum(x)})

# create a new assay to store ADT information
#NMF_assay <- CreateAssayObject(counts = t(usages))
NMF_assay <- CreateAssayObject(counts = usages_1)

# add this assay to the previously created Seurat object
DATA[["NMF"]] <- NMF_assay



DATA@active.assay = "NMF"

factors <- rownames(DATA)


allfactors = list()

colors = c("#FF0000","#00FF00","#0000FF")



for(i in factors){
  allfactors[[i]] = getGenePlots(i, topcolor = colors[as.numeric(gsub("X","",i))], qmax = 0.999, lowcolor = "#000000", imgalpha = 1, pt.alpha = 0.5)
}

plot_grid(plot_grid(plotlist =
lapply(1:length(allfactors[[1]]),function(j){
  plotcolors = lapply(allfactors, function(x){
    thisplot = x[[j]]
    ggplot_build(thisplot)$data[[1]]["fill"]
    })
  colordf = do.call("cbind",plotcolors)
  colors = apply(colordf, 1, function(crow){
    #do.call("rgb",as.list(apply(col2rgb(crow)^2,1, mean)/255^2))
    do.call("rgb",as.list(apply(col2rgb(crow)/255,1,max)))
    })
  g1 = allfactors[[1]][[j]]
  #g1[[x]]$data[,3] = as.factor(rownames(g1[[x]]$data))
  g1$data$cell = as.factor(rownames(g1$data))
  names(colors) = rownames(g1$data)
  g1$mapping$fill = g1$data$cell
  g1 = g1 + scale_fill_manual(values = colors) + NoLegend()
  return(g1)
  }), nrow = 1
),
plot_grid(get_legend(allfactors[[1]][[1]]),get_legend(allfactors[[2]][[1]]),
          get_legend(allfactors[[3]][[1]]), ncol = 3),
ncol = 1, rel_heights = c(3,1))

for(i in factors){
  inum = gsub("X","",i)
  topgenes = as.data.frame(t(sort(spectra[inum,], decreasing = TRUE)[1:15]))
  colnames(topgenes) = "Loading"
  topgenes$Gene = factor(rownames(topgenes), levels = rev(rownames(topgenes)))
  tgbar = ggplot(topgenes, aes(x = Loading, y = Gene)) + geom_bar(stat = "identity") + 
    theme(panel.background = element_blank(), panel.grid = element_blank()) + 
    labs(title = paste("Factor",inum))
  print(plot_grid(plotGene(i,  qmax = 0.999), 
                  plot_grid(tgbar,NULL, rel_heights = c(1.1,0.3), nrow = 2), 
                  rel_widths = c(3,1)))
}

DATA = SetIdent(DATA, value = "orig.ident")

d0 = data.frame(p.value = rep(NA,3), diff = rep(NA,3))
d3 = data.frame(p.value = rep(NA,3), diff = rep(NA,3))
for(i in 1:3){
  d0$p.value[i] = wilcox.test(DATA[["NMF"]]@counts[i,DATA$orig.ident=="GW_irrad0"],
              DATA[["NMF"]]@counts[i,DATA$orig.ident=="STD_irrad0"])$p.value
  d0$diff[i] = mean(DATA[["NMF"]]@counts[i,DATA$orig.ident=="GW_irrad0"]) -
              mean(DATA[["NMF"]]@counts[i,DATA$orig.ident=="STD_irrad0"])
  d3$p.value[i] = wilcox.test(DATA[["NMF"]]@counts[i,DATA$orig.ident=="GW_irrad3"],
              DATA[["NMF"]]@counts[i,DATA$orig.ident=="STD_irrad3"])$p.value
  d3$diff[i] = mean(DATA[["NMF"]]@counts[i,DATA$orig.ident=="GW_irrad3"]) -
              mean(DATA[["NMF"]]@counts[i,DATA$orig.ident=="STD_irrad3"])
}


plot_grid(plotlist = c(lapply(VlnPlot(DATA,features = factors, pt.size = 0, combine = FALSE), 
                            function(x){
                              x + NoLegend() + 
                                geom_point(alpha = 0.1, 
                                           position = position_jitter(width = 0.3),
                                           size = 0.1) +
                                geom_boxplot(width = 0.2, outlier.size = 0)}),
          list(DotPlot(DATA, features = c("X1","X2","X3")) + 
                 theme(legend.title = element_text(size = 8),
                       legend.text = element_text(size = 8))))
          )

GW vs STD d0

p-value Difference in mean
0e+00 -0.0508736
0e+00 0.0231163
2e-07 0.0277573

GW vs STD d3

p-value Difference in mean
0 -0.0065889
0 -0.0396660
0 0.0462549
topweights = as.data.frame(do.call("rbind",lapply(factors,function(i){
  inum = gsub("X","",i)
    
  df = as.data.frame(t(spectra[,names(sort(spectra[inum,], decreasing = TRUE)[1:20])]))
  df$gene = rownames(df)
  df

})))

topweights$gene = factor(topweights$gene, levels = topweights$gene)
heatmapdf = reshape2::melt(topweights)
colnames(heatmapdf) = c("Gene","Factor","Weight")
heatmapdf$Factor = factor(heatmapdf$Factor, levels = rev(c(1,2,3)))

ggplot(heatmapdf, aes(y = Factor, x = Gene)) +
  geom_tile(aes(fill = Weight)) + 
  theme(axis.text = element_text(size = 8),
        axis.title = element_text(size = 8),
        legend.text = element_text(size = 8),
        legend.title = element_text(size = 8),
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, size = 6),
        axis.ticks = element_blank(),
        panel.grid = element_blank(),
    legend.key.height = unit(10,"pt"),
    legend.key.width = unit(5,"pt"),
        panel.background = element_blank()) +
  scale_fill_viridis_c(option = "magma")

sessionInfo()
## R version 4.1.3 (2022-03-10)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: macOS Big Sur/Monterey 10.16
## 
## Matrix products: default
## BLAS/LAPACK: /Users/jenfra/miniconda3/envs/das_nature_2024_2/lib/libopenblasp-r0.3.27.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
## [1] cowplot_1.1.1       dplyr_1.1.2         ggplot2_3.4.2      
## [4] future.apply_1.11.0 future_1.32.0       Matrix_1.5-4.1     
## [7] SeuratObject_4.1.3  Seurat_4.3.0       
## 
## loaded via a namespace (and not attached):
##   [1] Rtsne_0.16             colorspace_2.1-0       deldir_1.0-9          
##   [4] ellipsis_0.3.2         ggridges_0.5.4         spatstat.data_3.0-1   
##   [7] farver_2.1.1           leiden_0.4.3           listenv_0.9.0         
##  [10] ggrepel_0.9.3          fansi_1.0.4            codetools_0.2-19      
##  [13] splines_4.1.3          cachem_1.0.8           knitr_1.43            
##  [16] polyclip_1.10-4        jsonlite_1.8.5         ica_1.0-3             
##  [19] cluster_2.1.4          png_0.1-8              uwot_0.1.14           
##  [22] shiny_1.7.4            sctransform_0.3.5      spatstat.sparse_3.0-1 
##  [25] compiler_4.1.3         httr_1.4.6             fastmap_1.1.1         
##  [28] lazyeval_0.2.2         cli_3.6.1              later_1.3.1           
##  [31] htmltools_0.5.5        tools_4.1.3            igraph_1.4.2          
##  [34] gtable_0.3.3           glue_1.6.2             RANN_2.6.1            
##  [37] reshape2_1.4.4         Rcpp_1.0.10            scattermore_1.1       
##  [40] jquerylib_0.1.4        vctrs_0.6.2            nlme_3.1-162          
##  [43] spatstat.explore_3.2-1 progressr_0.13.0       lmtest_0.9-40         
##  [46] spatstat.random_3.1-5  xfun_0.39              stringr_1.5.0         
##  [49] globals_0.16.2         mime_0.12              miniUI_0.1.1.1        
##  [52] lifecycle_1.0.3        irlba_2.3.5.1          goftest_1.2-3         
##  [55] MASS_7.3-58.3          zoo_1.8-12             scales_1.2.1          
##  [58] promises_1.2.0.1       spatstat.utils_3.0-3   parallel_4.1.3        
##  [61] RColorBrewer_1.1-3     yaml_2.3.7             reticulate_1.30       
##  [64] pbapply_1.7-0          gridExtra_2.3          sass_0.4.6            
##  [67] stringi_1.7.12         highr_0.10             rlang_1.1.1           
##  [70] pkgconfig_2.0.3        matrixStats_1.0.0      evaluate_0.21         
##  [73] lattice_0.21-8         ROCR_1.0-11            purrr_1.0.1           
##  [76] tensor_1.5             labeling_0.4.2         patchwork_1.1.2       
##  [79] htmlwidgets_1.6.2      tidyselect_1.2.0       parallelly_1.36.0     
##  [82] RcppAnnoy_0.0.20       plyr_1.8.8             magrittr_2.0.3        
##  [85] bookdown_0.34          R6_2.5.1               generics_0.1.3        
##  [88] withr_2.5.0            pillar_1.9.0           fitdistrplus_1.1-11   
##  [91] survival_3.5-5         abind_1.4-5            sp_1.6-1              
##  [94] tibble_3.2.1           KernSmooth_2.23-21     utf8_1.2.3            
##  [97] spatstat.geom_3.2-1    plotly_4.10.2          rmarkdown_2.22        
## [100] data.table_1.14.8      digest_0.6.31          xtable_1.8-4          
## [103] tidyr_1.3.0            httpuv_1.6.11          munsell_0.5.0         
## [106] viridisLite_0.4.1      bslib_0.5.0